#include<iostream>
#include<cstring>
#include "cmodel.h"
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/lexical_cast.hpp>

using boost::multiprecision::uint256_t;
using boost::multiprecision::uint128_t;
typedef unsigned char u8;
typedef unsigned int u32;

#define BYTELEN 8
#define BYTEM (char)(1 << BYTELEN)
#define BYTE1 (char)(1)
#define BIT(b, i) ((b >> (i)) & BYTE1)
#define FOURBITS(x) (x & 0xf)
#define BYTECOUNT(bitcount) ((bitcount + BYTELEN - 1) / BYTELEN)
#define OFFSET(ptr, unit_bitcount, offset) (ptr + BYTECOUNT(unit_bitcount) * offset)
#define MSB 0x80
#define LSB 0x01

// 将4个字节拼为32位
#define BUILD_UINT32(Byte3, Byte2, Byte1, Byte0)\
((u32)((u32)((Byte0) & 0x00FF)\
^ ((u32)((Byte1) & 0x00FF) << 8)\
^ ((u32)((Byte2) & 0x00FF) << 16)\
^ ((u32)((Byte3) & 0x00FF) << 24)))
#define PUTU32(ct, st) { (ct)[0] = (u8)((st) >> 24); (ct)[1] = (u8)((st) >> 16); (ct)[2] = (u8)((st) >>  8); (ct)[3] = (u8)(st);}

#define MON_MUL_BLOCK_LENGTH 128  // 硬件实现大整数乘法时是分块进行的，此变量为每块的位宽
#define ADD_BLOCK_LENGTH 128  // 加法操作的每个单元的位宽
#define AES_BLOCK_SIZE 16
static const __uint128_t MAX_128 = -1;
static const uint256_t MAX_256 = -1;

static const u32 Te0[256] = {
        0xc66363a5, 0xf87c7c84, 0xee777799, 0xf67b7b8d, 0xfff2f20d, 0xd66b6bbd, 0xde6f6fb1, 0x91c5c554,
        0x60303050, 0x02010103, 0xce6767a9, 0x562b2b7d, 0xe7fefe19, 0xb5d7d762, 0x4dababe6, 0xec76769a,

        0x8fcaca45, 0x1f82829d, 0x89c9c940, 0xfa7d7d87, 0xeffafa15, 0xb25959eb, 0x8e4747c9, 0xfbf0f00b,
        0x41adadec, 0xb3d4d467, 0x5fa2a2fd, 0x45afafea, 0x239c9cbf, 0x53a4a4f7, 0xe4727296, 0x9bc0c05b,

        0x75b7b7c2, 0xe1fdfd1c, 0x3d9393ae, 0x4c26266a, 0x6c36365a, 0x7e3f3f41, 0xf5f7f702, 0x83cccc4f,
        0x6834345c, 0x51a5a5f4, 0xd1e5e534, 0xf9f1f108, 0xe2717193, 0xabd8d873, 0x62313153, 0x2a15153f,

        0x0804040c, 0x95c7c752, 0x46232365, 0x9dc3c35e, 0x30181828, 0x379696a1, 0x0a05050f, 0x2f9a9ab5,
        0x0e070709, 0x24121236, 0x1b80809b, 0xdfe2e23d, 0xcdebeb26, 0x4e272769, 0x7fb2b2cd, 0xea75759f,

        0x1209091b, 0x1d83839e, 0x582c2c74, 0x341a1a2e, 0x361b1b2d, 0xdc6e6eb2, 0xb45a5aee, 0x5ba0a0fb,
        0xa45252f6, 0x763b3b4d, 0xb7d6d661, 0x7db3b3ce, 0x5229297b, 0xdde3e33e, 0x5e2f2f71, 0x13848497,

        0xa65353f5, 0xb9d1d168, 0x00000000, 0xc1eded2c, 0x40202060, 0xe3fcfc1f, 0x79b1b1c8, 0xb65b5bed,
        0xd46a6abe, 0x8dcbcb46, 0x67bebed9, 0x7239394b, 0x944a4ade, 0x984c4cd4, 0xb05858e8, 0x85cfcf4a,

        0xbbd0d06b, 0xc5efef2a, 0x4faaaae5, 0xedfbfb16, 0x864343c5, 0x9a4d4dd7, 0x66333355, 0x11858594,
        0x8a4545cf, 0xe9f9f910, 0x04020206, 0xfe7f7f81, 0xa05050f0, 0x783c3c44, 0x259f9fba, 0x4ba8a8e3,

        0xa25151f3, 0x5da3a3fe, 0x804040c0, 0x058f8f8a, 0x3f9292ad, 0x219d9dbc, 0x70383848, 0xf1f5f504,
        0x63bcbcdf, 0x77b6b6c1, 0xafdada75, 0x42212163, 0x20101030, 0xe5ffff1a, 0xfdf3f30e, 0xbfd2d26d,

        0x81cdcd4c, 0x180c0c14, 0x26131335, 0xc3ecec2f, 0xbe5f5fe1, 0x359797a2, 0x884444cc, 0x2e171739,
        0x93c4c457, 0x55a7a7f2, 0xfc7e7e82, 0x7a3d3d47, 0xc86464ac, 0xba5d5de7, 0x3219192b, 0xe6737395,

        0xc06060a0, 0x19818198, 0x9e4f4fd1, 0xa3dcdc7f, 0x44222266, 0x542a2a7e, 0x3b9090ab, 0x0b888883,
        0x8c4646ca, 0xc7eeee29, 0x6bb8b8d3, 0x2814143c, 0xa7dede79, 0xbc5e5ee2, 0x160b0b1d, 0xaddbdb76,

        0xdbe0e03b, 0x64323256, 0x743a3a4e, 0x140a0a1e, 0x924949db, 0x0c06060a, 0x4824246c, 0xb85c5ce4,
        0x9fc2c25d, 0xbdd3d36e, 0x43acacef, 0xc46262a6, 0x399191a8, 0x319595a4, 0xd3e4e437, 0xf279798b,

        0xd5e7e732, 0x8bc8c843, 0x6e373759, 0xda6d6db7, 0x018d8d8c, 0xb1d5d564, 0x9c4e4ed2, 0x49a9a9e0,
        0xd86c6cb4, 0xac5656fa, 0xf3f4f407, 0xcfeaea25, 0xca6565af, 0xf47a7a8e, 0x47aeaee9, 0x10080818,

        0x6fbabad5, 0xf0787888, 0x4a25256f, 0x5c2e2e72, 0x381c1c24, 0x57a6a6f1, 0x73b4b4c7, 0x97c6c651,
        0xcbe8e823, 0xa1dddd7c, 0xe874749c, 0x3e1f1f21, 0x964b4bdd, 0x61bdbddc, 0x0d8b8b86, 0x0f8a8a85,

        0xe0707090, 0x7c3e3e42, 0x71b5b5c4, 0xcc6666aa, 0x904848d8, 0x06030305, 0xf7f6f601, 0x1c0e0e12,
        0xc26161a3, 0x6a35355f, 0xae5757f9, 0x69b9b9d0, 0x17868691, 0x99c1c158, 0x3a1d1d27, 0x279e9eb9,

        0xd9e1e138, 0xebf8f813, 0x2b9898b3, 0x22111133, 0xd26969bb, 0xa9d9d970, 0x078e8e89, 0x339494a7,
        0x2d9b9bb6, 0x3c1e1e22, 0x15878792, 0xc9e9e920, 0x87cece49, 0xaa5555ff, 0x50282878, 0xa5dfdf7a,

        0x038c8c8f, 0x59a1a1f8, 0x09898980, 0x1a0d0d17, 0x65bfbfda, 0xd7e6e631, 0x844242c6, 0xd06868b8,
        0x824141c3, 0x299999b0, 0x5a2d2d77, 0x1e0f0f11, 0x7bb0b0cb, 0xa85454fc, 0x6dbbbbd6, 0x2c16163a,
};

static const u32 Te1[256] = {
        0xa5c66363, 0x84f87c7c, 0x99ee7777, 0x8df67b7b, 0x0dfff2f2, 0xbdd66b6b, 0xb1de6f6f, 0x5491c5c5,
        0x50603030, 0x03020101, 0xa9ce6767, 0x7d562b2b, 0x19e7fefe, 0x62b5d7d7, 0xe64dabab, 0x9aec7676,

        0x458fcaca, 0x9d1f8282, 0x4089c9c9, 0x87fa7d7d, 0x15effafa, 0xebb25959, 0xc98e4747, 0x0bfbf0f0,
        0xec41adad, 0x67b3d4d4, 0xfd5fa2a2, 0xea45afaf, 0xbf239c9c, 0xf753a4a4, 0x96e47272, 0x5b9bc0c0,

        0xc275b7b7, 0x1ce1fdfd, 0xae3d9393, 0x6a4c2626, 0x5a6c3636, 0x417e3f3f, 0x02f5f7f7, 0x4f83cccc,
        0x5c683434, 0xf451a5a5, 0x34d1e5e5, 0x08f9f1f1, 0x93e27171, 0x73abd8d8, 0x53623131, 0x3f2a1515,

        0x0c080404, 0x5295c7c7, 0x65462323, 0x5e9dc3c3, 0x28301818, 0xa1379696, 0x0f0a0505, 0xb52f9a9a,
        0x090e0707, 0x36241212, 0x9b1b8080, 0x3ddfe2e2, 0x26cdebeb, 0x694e2727, 0xcd7fb2b2, 0x9fea7575,

        0x1b120909, 0x9e1d8383, 0x74582c2c, 0x2e341a1a, 0x2d361b1b, 0xb2dc6e6e, 0xeeb45a5a, 0xfb5ba0a0,
        0xf6a45252, 0x4d763b3b, 0x61b7d6d6, 0xce7db3b3, 0x7b522929, 0x3edde3e3, 0x715e2f2f, 0x97138484,

        0xf5a65353, 0x68b9d1d1, 0x00000000, 0x2cc1eded, 0x60402020, 0x1fe3fcfc, 0xc879b1b1, 0xedb65b5b,
        0xbed46a6a, 0x468dcbcb, 0xd967bebe, 0x4b723939, 0xde944a4a, 0xd4984c4c, 0xe8b05858, 0x4a85cfcf,

        0x6bbbd0d0, 0x2ac5efef, 0xe54faaaa, 0x16edfbfb, 0xc5864343, 0xd79a4d4d, 0x55663333, 0x94118585,
        0xcf8a4545, 0x10e9f9f9, 0x06040202, 0x81fe7f7f, 0xf0a05050, 0x44783c3c, 0xba259f9f, 0xe34ba8a8,

        0xf3a25151, 0xfe5da3a3, 0xc0804040, 0x8a058f8f, 0xad3f9292, 0xbc219d9d, 0x48703838, 0x04f1f5f5,
        0xdf63bcbc, 0xc177b6b6, 0x75afdada, 0x63422121, 0x30201010, 0x1ae5ffff, 0x0efdf3f3, 0x6dbfd2d2,

        0x4c81cdcd, 0x14180c0c, 0x35261313, 0x2fc3ecec, 0xe1be5f5f, 0xa2359797, 0xcc884444, 0x392e1717,
        0x5793c4c4, 0xf255a7a7, 0x82fc7e7e, 0x477a3d3d, 0xacc86464, 0xe7ba5d5d, 0x2b321919, 0x95e67373,

        0xa0c06060, 0x98198181, 0xd19e4f4f, 0x7fa3dcdc, 0x66442222, 0x7e542a2a, 0xab3b9090, 0x830b8888,
        0xca8c4646, 0x29c7eeee, 0xd36bb8b8, 0x3c281414, 0x79a7dede, 0xe2bc5e5e, 0x1d160b0b, 0x76addbdb,

        0x3bdbe0e0, 0x56643232, 0x4e743a3a, 0x1e140a0a, 0xdb924949, 0x0a0c0606, 0x6c482424, 0xe4b85c5c,
        0x5d9fc2c2, 0x6ebdd3d3, 0xef43acac, 0xa6c46262, 0xa8399191, 0xa4319595, 0x37d3e4e4, 0x8bf27979,

        0x32d5e7e7, 0x438bc8c8, 0x596e3737, 0xb7da6d6d, 0x8c018d8d, 0x64b1d5d5, 0xd29c4e4e, 0xe049a9a9,
        0xb4d86c6c, 0xfaac5656, 0x07f3f4f4, 0x25cfeaea, 0xafca6565, 0x8ef47a7a, 0xe947aeae, 0x18100808,

        0xd56fbaba, 0x88f07878, 0x6f4a2525, 0x725c2e2e, 0x24381c1c, 0xf157a6a6, 0xc773b4b4, 0x5197c6c6,
        0x23cbe8e8, 0x7ca1dddd, 0x9ce87474, 0x213e1f1f, 0xdd964b4b, 0xdc61bdbd, 0x860d8b8b, 0x850f8a8a,

        0x90e07070, 0x427c3e3e, 0xc471b5b5, 0xaacc6666, 0xd8904848, 0x05060303, 0x01f7f6f6, 0x121c0e0e,
        0xa3c26161, 0x5f6a3535, 0xf9ae5757, 0xd069b9b9, 0x91178686, 0x5899c1c1, 0x273a1d1d, 0xb9279e9e,

        0x38d9e1e1, 0x13ebf8f8, 0xb32b9898, 0x33221111, 0xbbd26969, 0x70a9d9d9, 0x89078e8e, 0xa7339494,
        0xb62d9b9b, 0x223c1e1e, 0x92158787, 0x20c9e9e9, 0x4987cece, 0xffaa5555, 0x78502828, 0x7aa5dfdf,

        0x8f038c8c, 0xf859a1a1, 0x80098989, 0x171a0d0d, 0xda65bfbf, 0x31d7e6e6, 0xc6844242, 0xb8d06868,
        0xc3824141, 0xb0299999, 0x775a2d2d, 0x111e0f0f, 0xcb7bb0b0, 0xfca85454, 0xd66dbbbb, 0x3a2c1616,
};

static const u32 Te2[256] = {
        0x63a5c663, 0x7c84f87c, 0x7799ee77, 0x7b8df67b, 0xf20dfff2, 0x6bbdd66b, 0x6fb1de6f, 0xc55491c5,
        0x30506030, 0x01030201, 0x67a9ce67, 0x2b7d562b, 0xfe19e7fe, 0xd762b5d7, 0xabe64dab, 0x769aec76,

        0xca458fca, 0x829d1f82, 0xc94089c9, 0x7d87fa7d, 0xfa15effa, 0x59ebb259, 0x47c98e47, 0xf00bfbf0,
        0xadec41ad, 0xd467b3d4, 0xa2fd5fa2, 0xafea45af, 0x9cbf239c, 0xa4f753a4, 0x7296e472, 0xc05b9bc0,

        0xb7c275b7, 0xfd1ce1fd, 0x93ae3d93, 0x266a4c26, 0x365a6c36, 0x3f417e3f, 0xf702f5f7, 0xcc4f83cc,
        0x345c6834, 0xa5f451a5, 0xe534d1e5, 0xf108f9f1, 0x7193e271, 0xd873abd8, 0x31536231, 0x153f2a15,

        0x040c0804, 0xc75295c7, 0x23654623, 0xc35e9dc3, 0x18283018, 0x96a13796, 0x050f0a05, 0x9ab52f9a,
        0x07090e07, 0x12362412, 0x809b1b80, 0xe23ddfe2, 0xeb26cdeb, 0x27694e27, 0xb2cd7fb2, 0x759fea75,

        0x091b1209, 0x839e1d83, 0x2c74582c, 0x1a2e341a, 0x1b2d361b, 0x6eb2dc6e, 0x5aeeb45a, 0xa0fb5ba0,
        0x52f6a452, 0x3b4d763b, 0xd661b7d6, 0xb3ce7db3, 0x297b5229, 0xe33edde3, 0x2f715e2f, 0x84971384,

        0x53f5a653, 0xd168b9d1, 0x00000000, 0xed2cc1ed, 0x20604020, 0xfc1fe3fc, 0xb1c879b1, 0x5bedb65b,
        0x6abed46a, 0xcb468dcb, 0xbed967be, 0x394b7239, 0x4ade944a, 0x4cd4984c, 0x58e8b058, 0xcf4a85cf,

        0xd06bbbd0, 0xef2ac5ef, 0xaae54faa, 0xfb16edfb, 0x43c58643, 0x4dd79a4d, 0x33556633, 0x85941185,
        0x45cf8a45, 0xf910e9f9, 0x02060402, 0x7f81fe7f, 0x50f0a050, 0x3c44783c, 0x9fba259f, 0xa8e34ba8,

        0x51f3a251, 0xa3fe5da3, 0x40c08040, 0x8f8a058f, 0x92ad3f92, 0x9dbc219d, 0x38487038, 0xf504f1f5,
        0xbcdf63bc, 0xb6c177b6, 0xda75afda, 0x21634221, 0x10302010, 0xff1ae5ff, 0xf30efdf3, 0xd26dbfd2,

        0xcd4c81cd, 0x0c14180c, 0x13352613, 0xec2fc3ec, 0x5fe1be5f, 0x97a23597, 0x44cc8844, 0x17392e17,
        0xc45793c4, 0xa7f255a7, 0x7e82fc7e, 0x3d477a3d, 0x64acc864, 0x5de7ba5d, 0x192b3219, 0x7395e673,

        0x60a0c060, 0x81981981, 0x4fd19e4f, 0xdc7fa3dc, 0x22664422, 0x2a7e542a, 0x90ab3b90, 0x88830b88,
        0x46ca8c46, 0xee29c7ee, 0xb8d36bb8, 0x143c2814, 0xde79a7de, 0x5ee2bc5e, 0x0b1d160b, 0xdb76addb,

        0xe03bdbe0, 0x32566432, 0x3a4e743a, 0x0a1e140a, 0x49db9249, 0x060a0c06, 0x246c4824, 0x5ce4b85c,
        0xc25d9fc2, 0xd36ebdd3, 0xacef43ac, 0x62a6c462, 0x91a83991, 0x95a43195, 0xe437d3e4, 0x798bf279,

        0xe732d5e7, 0xc8438bc8, 0x37596e37, 0x6db7da6d, 0x8d8c018d, 0xd564b1d5, 0x4ed29c4e, 0xa9e049a9,
        0x6cb4d86c, 0x56faac56, 0xf407f3f4, 0xea25cfea, 0x65afca65, 0x7a8ef47a, 0xaee947ae, 0x08181008,

        0xbad56fba, 0x7888f078, 0x256f4a25, 0x2e725c2e, 0x1c24381c, 0xa6f157a6, 0xb4c773b4, 0xc65197c6,
        0xe823cbe8, 0xdd7ca1dd, 0x749ce874, 0x1f213e1f, 0x4bdd964b, 0xbddc61bd, 0x8b860d8b, 0x8a850f8a,

        0x7090e070, 0x3e427c3e, 0xb5c471b5, 0x66aacc66, 0x48d89048, 0x03050603, 0xf601f7f6, 0x0e121c0e,
        0x61a3c261, 0x355f6a35, 0x57f9ae57, 0xb9d069b9, 0x86911786, 0xc15899c1, 0x1d273a1d, 0x9eb9279e,

        0xe138d9e1, 0xf813ebf8, 0x98b32b98, 0x11332211, 0x69bbd269, 0xd970a9d9, 0x8e89078e, 0x94a73394,
        0x9bb62d9b, 0x1e223c1e, 0x87921587, 0xe920c9e9, 0xce4987ce, 0x55ffaa55, 0x28785028, 0xdf7aa5df,

        0x8c8f038c, 0xa1f859a1, 0x89800989, 0x0d171a0d, 0xbfda65bf, 0xe631d7e6, 0x42c68442, 0x68b8d068,
        0x41c38241, 0x99b02999, 0x2d775a2d, 0x0f111e0f, 0xb0cb7bb0, 0x54fca854, 0xbbd66dbb, 0x163a2c16,
};

static const u32 Te3[256] = {
        0x6363a5c6, 0x7c7c84f8, 0x777799ee, 0x7b7b8df6, 0xf2f20dff, 0x6b6bbdd6, 0x6f6fb1de, 0xc5c55491,
        0x30305060, 0x01010302, 0x6767a9ce, 0x2b2b7d56, 0xfefe19e7, 0xd7d762b5, 0xababe64d, 0x76769aec,

        0xcaca458f, 0x82829d1f, 0xc9c94089, 0x7d7d87fa, 0xfafa15ef, 0x5959ebb2, 0x4747c98e, 0xf0f00bfb,
        0xadadec41, 0xd4d467b3, 0xa2a2fd5f, 0xafafea45, 0x9c9cbf23, 0xa4a4f753, 0x727296e4, 0xc0c05b9b,

        0xb7b7c275, 0xfdfd1ce1, 0x9393ae3d, 0x26266a4c, 0x36365a6c, 0x3f3f417e, 0xf7f702f5, 0xcccc4f83,
        0x34345c68, 0xa5a5f451, 0xe5e534d1, 0xf1f108f9, 0x717193e2, 0xd8d873ab, 0x31315362, 0x15153f2a,

        0x04040c08, 0xc7c75295, 0x23236546, 0xc3c35e9d, 0x18182830, 0x9696a137, 0x05050f0a, 0x9a9ab52f,
        0x0707090e, 0x12123624, 0x80809b1b, 0xe2e23ddf, 0xebeb26cd, 0x2727694e, 0xb2b2cd7f, 0x75759fea,

        0x09091b12, 0x83839e1d, 0x2c2c7458, 0x1a1a2e34, 0x1b1b2d36, 0x6e6eb2dc, 0x5a5aeeb4, 0xa0a0fb5b,
        0x5252f6a4, 0x3b3b4d76, 0xd6d661b7, 0xb3b3ce7d, 0x29297b52, 0xe3e33edd, 0x2f2f715e, 0x84849713,

        0x5353f5a6, 0xd1d168b9, 0x00000000, 0xeded2cc1, 0x20206040, 0xfcfc1fe3, 0xb1b1c879, 0x5b5bedb6,
        0x6a6abed4, 0xcbcb468d, 0xbebed967, 0x39394b72, 0x4a4ade94, 0x4c4cd498, 0x5858e8b0, 0xcfcf4a85,

        0xd0d06bbb, 0xefef2ac5, 0xaaaae54f, 0xfbfb16ed, 0x4343c586, 0x4d4dd79a, 0x33335566, 0x85859411,
        0x4545cf8a, 0xf9f910e9, 0x02020604, 0x7f7f81fe, 0x5050f0a0, 0x3c3c4478, 0x9f9fba25, 0xa8a8e34b,

        0x5151f3a2, 0xa3a3fe5d, 0x4040c080, 0x8f8f8a05, 0x9292ad3f, 0x9d9dbc21, 0x38384870, 0xf5f504f1,
        0xbcbcdf63, 0xb6b6c177, 0xdada75af, 0x21216342, 0x10103020, 0xffff1ae5, 0xf3f30efd, 0xd2d26dbf,

        0xcdcd4c81, 0x0c0c1418, 0x13133526, 0xecec2fc3, 0x5f5fe1be, 0x9797a235, 0x4444cc88, 0x1717392e,
        0xc4c45793, 0xa7a7f255, 0x7e7e82fc, 0x3d3d477a, 0x6464acc8, 0x5d5de7ba, 0x19192b32, 0x737395e6,

        0x6060a0c0, 0x81819819, 0x4f4fd19e, 0xdcdc7fa3, 0x22226644, 0x2a2a7e54, 0x9090ab3b, 0x8888830b,
        0x4646ca8c, 0xeeee29c7, 0xb8b8d36b, 0x14143c28, 0xdede79a7, 0x5e5ee2bc, 0x0b0b1d16, 0xdbdb76ad,

        0xe0e03bdb, 0x32325664, 0x3a3a4e74, 0x0a0a1e14, 0x4949db92, 0x06060a0c, 0x24246c48, 0x5c5ce4b8,
        0xc2c25d9f, 0xd3d36ebd, 0xacacef43, 0x6262a6c4, 0x9191a839, 0x9595a431, 0xe4e437d3, 0x79798bf2,

        0xe7e732d5, 0xc8c8438b, 0x3737596e, 0x6d6db7da, 0x8d8d8c01, 0xd5d564b1, 0x4e4ed29c, 0xa9a9e049,
        0x6c6cb4d8, 0x5656faac, 0xf4f407f3, 0xeaea25cf, 0x6565afca, 0x7a7a8ef4, 0xaeaee947, 0x08081810,

        0xbabad56f, 0x787888f0, 0x25256f4a, 0x2e2e725c, 0x1c1c2438, 0xa6a6f157, 0xb4b4c773, 0xc6c65197,
        0xe8e823cb, 0xdddd7ca1, 0x74749ce8, 0x1f1f213e, 0x4b4bdd96, 0xbdbddc61, 0x8b8b860d, 0x8a8a850f,

        0x707090e0, 0x3e3e427c, 0xb5b5c471, 0x6666aacc, 0x4848d890, 0x03030506, 0xf6f601f7, 0x0e0e121c,
        0x6161a3c2, 0x35355f6a, 0x5757f9ae, 0xb9b9d069, 0x86869117, 0xc1c15899, 0x1d1d273a, 0x9e9eb927,

        0xe1e138d9, 0xf8f813eb, 0x9898b32b, 0x11113322, 0x6969bbd2, 0xd9d970a9, 0x8e8e8907, 0x9494a733,
        0x9b9bb62d, 0x1e1e223c, 0x87879215, 0xe9e920c9, 0xcece4987, 0x5555ffaa, 0x28287850, 0xdfdf7aa5,

        0x8c8c8f03, 0xa1a1f859, 0x89898009, 0x0d0d171a, 0xbfbfda65, 0xe6e631d7, 0x4242c684, 0x6868b8d0,
        0x4141c382, 0x9999b029, 0x2d2d775a, 0x0f0f111e, 0xb0b0cb7b, 0x5454fca8, 0xbbbbd66d, 0x16163a2c,
};

static const u8 Sbox[256] = {
        0x63, 0x7c, 0x77, 0x7b, 0xf2, 0x6b, 0x6f, 0xc5, 0x30, 0x01, 0x67, 0x2b, 0xfe, 0xd7, 0xab, 0x76,
        0xca, 0x82, 0xc9, 0x7d, 0xfa, 0x59, 0x47, 0xf0, 0xad, 0xd4, 0xa2, 0xaf, 0x9c, 0xa4, 0x72, 0xc0,
        0xb7, 0xfd, 0x93, 0x26, 0x36, 0x3f, 0xf7, 0xcc, 0x34, 0xa5, 0xe5, 0xf1, 0x71, 0xd8, 0x31, 0x15,
        0x04, 0xc7, 0x23, 0xc3, 0x18, 0x96, 0x05, 0x9a, 0x07, 0x12, 0x80, 0xe2, 0xeb, 0x27, 0xb2, 0x75,
        0x09, 0x83, 0x2c, 0x1a, 0x1b, 0x6e, 0x5a, 0xa0, 0x52, 0x3b, 0xd6, 0xb3, 0x29, 0xe3, 0x2f, 0x84,
        0x53, 0xd1, 0x00, 0xed, 0x20, 0xfc, 0xb1, 0x5b, 0x6a, 0xcb, 0xbe, 0x39, 0x4a, 0x4c, 0x58, 0xcf,
        0xd0, 0xef, 0xaa, 0xfb, 0x43, 0x4d, 0x33, 0x85, 0x45, 0xf9, 0x02, 0x7f, 0x50, 0x3c, 0x9f, 0xa8,
        0x51, 0xa3, 0x40, 0x8f, 0x92, 0x9d, 0x38, 0xf5, 0xbc, 0xb6, 0xda, 0x21, 0x10, 0xff, 0xf3, 0xd2,
        0xcd, 0x0c, 0x13, 0xec, 0x5f, 0x97, 0x44, 0x17, 0xc4, 0xa7, 0x7e, 0x3d, 0x64, 0x5d, 0x19, 0x73,
        0x60, 0x81, 0x4f, 0xdc, 0x22, 0x2a, 0x90, 0x88, 0x46, 0xee, 0xb8, 0x14, 0xde, 0x5e, 0x0b, 0xdb,
        0xe0, 0x32, 0x3a, 0x0a, 0x49, 0x06, 0x24, 0x5c, 0xc2, 0xd3, 0xac, 0x62, 0x91, 0x95, 0xe4, 0x79,
        0xe7, 0xc8, 0x37, 0x6d, 0x8d, 0xd5, 0x4e, 0xa9, 0x6c, 0x56, 0xf4, 0xea, 0x65, 0x7a, 0xae, 0x08,
        0xba, 0x78, 0x25, 0x2e, 0x1c, 0xa6, 0xb4, 0xc6, 0xe8, 0xdd, 0x74, 0x1f, 0x4b, 0xbd, 0x8b, 0x8a,
        0x70, 0x3e, 0xb5, 0x66, 0x48, 0x03, 0xf6, 0x0e, 0x61, 0x35, 0x57, 0xb9, 0x86, 0xc1, 0x1d, 0x9e,
        0xe1, 0xf8, 0x98, 0x11, 0x69, 0xd9, 0x8e, 0x94, 0x9b, 0x1e, 0x87, 0xe9, 0xce, 0x55, 0x28, 0xdf,
        0x8c, 0xa1, 0x89, 0x0d, 0xbf, 0xe6, 0x42, 0x68, 0x41, 0x99, 0x2d, 0x0f, 0xb0, 0x54, 0xbb, 0x16 };

static const u32 rcon[] = {
        0x01000000, 0x02000000, 0x04000000, 0x08000000,
        0x10000000, 0x20000000, 0x40000000, 0x80000000,
        0x1B000000, 0x36000000,
};

uint16_t checksum(unsigned char *data, size_t byte_len) {
    uint32_t cksum = 0;
    /* use u8_data_ptr and u16_data_ptr*/
    //uint16_t temp;
    // const uint8_t *u8_data_ptr = reinterpret_cast<const uint8_t*>(data);
    // while (byte_len > 1) {
    //   temp = static_cast<uint16_t>(*u8_data_ptr++) << 8;
    //   temp |= (*u8_data_ptr++);
    //   cksum += temp;
    //   byte_len -= 2;
    // }
    // if (byte_len) {
    //   cksum += *(u8_data_ptr) << 8;
    // }
    const uint16_t *u16_data_ptr = reinterpret_cast<const uint16_t*>(data);
    while (byte_len > 1) {
        cksum += *u16_data_ptr++;
        byte_len -= 2;
    }
    // if byte_len is odd
    if (byte_len) {
        cksum += *(reinterpret_cast<const uint8_t*>(u16_data_ptr));
    }
    // 将32位cksum转换为16位, 若高16bit位有进位，进位到低16bit位
    // 若将cksum表示为0xabcdefgh, 这步代表0x0000abcd + 0x0000efgh
    cksum = (cksum >> 16) + (cksum & 0xffff);
    // 若仍有进位(且最多1位)得到0x0001xyzt, 这步继续将进位相加 0x000000001 + 0x0001xyzt
    cksum += (cksum >> 16);
    // 0x000000001 + 0x0001xyzt不可能再次进位到高16bit
    //取反, uint_16t转换将抛弃高16bit位
    return static_cast<uint16_t>(~cksum);
}

void _encode_int_to_char(const void *in, size_t len, char *out) {
    const char *tmpin = (const char *)in;
    for (size_t i = 0; i < len; i++) {
        out[len - 1 - i] = tmpin[i];
    }
}

void _decode_char_to_int(const char *in, size_t len, void *out) {
    char *tmpout = (char *)out;
    for (size_t i = 0; i < len; i++) {
        tmpout[i] = in[len - 1 - i];
    }
}

void _encode_uint256_t_to_char(const uint256_t *in, char *out) {
    uint128_t high = boost::lexical_cast<uint128_t>(*in >> 128);
    uint128_t low = boost::lexical_cast<uint128_t>(*in);
    char *high_char = (char*)(&high);
    char *low_char = (char*)(&low);
    for (size_t i = 0; i < 16; i++) {
        out[15 - i] = high_char[i];
        out[31 - i] = low_char[i];
    }
}

void _decode_char_to_uint256_t(const char *in, uint256_t *out) {
    uint128_t high;
    uint128_t low;
    char *high_char = (char*)(&high);
    char *low_char = (char*)(&low);
    for (size_t i = 0; i < 16; i++) {
        high_char[i] = in[15 - i];
        low_char[i] = in[31 - i];
    }
    uint256_t high_256 =(boost::lexical_cast<uint256_t>(high)) << 128;
    uint256_t low_256 =boost::lexical_cast<uint256_t>(low);
    *out = high_256 + low_256;
}

void _multi_64(const char *a, const char *b, char *result) {
    uint64_t a_temp = 0;
    uint64_t b_temp = 0;
    _decode_char_to_int(a, sizeof(a_temp), &a_temp);
    _decode_char_to_int(b, sizeof(b_temp), &b_temp);
    __uint128_t result_temp = 0;
    result_temp = (__uint128_t)(a_temp) * (__uint128_t)(b_temp);
    _encode_int_to_char(&result_temp, sizeof(__uint128_t), result);
}

void _multi_128(const char *a, const char *b, char *result) {
    __uint128_t a_temp = 0;
    __uint128_t b_temp = 0;
    _decode_char_to_int(a, 16, &a_temp);
    _decode_char_to_int(b, 16, &b_temp);
    uint256_t result_temp = 0;
    result_temp = (uint256_t)(a_temp) * (uint256_t)(b_temp);
    _encode_uint256_t_to_char(&result_temp, result);
}

void _add_128(const char *a, const char *b, const int carry_in, char *result, bool *carry_out) {
    __uint128_t a_temp = 0;
    __uint128_t b_temp = 0;
    _decode_char_to_int(a, sizeof(a_temp), &a_temp);
    _decode_char_to_int(b, sizeof(b_temp), &b_temp);
    *carry_out = a_temp > MAX_128 - carry_in;
    a_temp = a_temp + carry_in;
    *carry_out = *carry_out | (a_temp > MAX_128 - b_temp);  // if a_temp + b_temp + carry_in > MAX_128, then set carry=true
    a_temp = a_temp + b_temp;
    _encode_int_to_char(&a_temp, sizeof(a_temp), result);
}

void _add_256(const char *a, const char *b, const int carry_in, char *result, bool *carry_out) {
    uint256_t a_temp = 0;
    uint256_t b_temp = 0;
    _decode_char_to_uint256_t(a, &a_temp);
    _decode_char_to_uint256_t(b, &b_temp);
    *carry_out = a_temp > MAX_256 - carry_in;
    a_temp = a_temp + carry_in;
    *carry_out = *carry_out | (a_temp > MAX_256 - b_temp);  // if a_temp + b_temp + carry_in > MAX_128, then set carry=true
    a_temp = a_temp + b_temp;
    _encode_uint256_t_to_char(&a_temp, result);
}

void _sub_128(const char *a, const char *b, const int carry_in, char *result, bool *carry_out) {
    __uint128_t a_temp = 0;
    __uint128_t b_temp = 0;
    _decode_char_to_int(a, sizeof(a_temp), &a_temp);
    _decode_char_to_int(b, sizeof(b_temp), &b_temp);
    *carry_out = a_temp < (unsigned int)carry_in;
    a_temp = a_temp - carry_in;
    *carry_out = *carry_out | (a_temp < b_temp);  // if a_temp - carry_in < b_temp, then set carry=true
    a_temp = a_temp - b_temp;
    _encode_int_to_char(&a_temp, sizeof(a_temp), result);
}

/*
 * @brief: 模加单元，目标是输出 a+b mod p，其中0 <= a,b < p
 * @param[in][out] a: 加数，同时也是存放的输出结果
 * @param[in] b: 加数
 * @param[in] p: 模数
 * @param[in] p_bitcount: 模数的比特数，确保是8的倍数，一般而言是2048或4096
 */
void _mod_add_unit(char *a, char *b, char *p, int p_bitcount) {
    if (p_bitcount < ADD_BLOCK_LENGTH) {
        __uint128_t a_temp = 0;
        __uint128_t b_temp = 0;
        _decode_char_to_int(a, BYTECOUNT(p_bitcount), &a_temp);
        _decode_char_to_int(b, BYTECOUNT(p_bitcount), &b_temp);
        a_temp += b_temp;
        _decode_char_to_int(p, BYTECOUNT(p_bitcount), &b_temp);
        if (a_temp > b_temp) {
            a_temp = a_temp - b_temp;  // if a + b > p, then the result is a + b - p
        }
        _encode_int_to_char(&a_temp, BYTECOUNT(p_bitcount), a);
        return;
    }

    // ADD_BLOCK_LENGTH = 128, p_bitcount = 2048, block_num = 16
    int block_count = (p_bitcount + ADD_BLOCK_LENGTH - 1) / ADD_BLOCK_LENGTH;
    int block_size = BYTECOUNT(ADD_BLOCK_LENGTH);
    // 在加减法实现中，carry_in和carry_out不能指向同一块空间
    int carry_in = 0;
    bool carry_out = false;
    for (int i = 0; i < block_count; i++) {
        _add_128(a + (block_count - 1 - i) * block_size, b + (block_count - 1 - i) * block_size, carry_in,
                 a + (block_count - 1 - i) * block_size, &carry_out);
        // process the carry
        carry_in = (int)carry_out;
    }

    // compare the size of result and p
    bool flag = false;
    if (carry_out) {
        flag = true;
    }
    if (!flag) {
        for (int i = 0; i < block_size * block_count; i++) {
            if ((int)((unsigned char)(a[i])) > (int)((unsigned char)(p[i]))) {
                flag = true;
                break;
            }
            else if ((int)((unsigned char)(a[i])) < (int)((unsigned char)(p[i]))) {
                flag = false;
                break;
            }
            else {
                continue;
            }
        }
    }

    // if result > p, then compute output = result - p
    carry_in = 0;
    if (flag) {
        for (int i = 0; i < block_count; i++) {
            _sub_128(a + (block_count - 1 - i) * block_size, p + (block_count - 1 - i) * block_size, carry_in,
                     a + (block_count - 1 - i) * block_size, &carry_out);
            // process the carry
            carry_in = (int)carry_out;
        }
    }
}

/*
 * @brief: 蒙哥马利模乘单元，目标是输出 a*b*r^(-1) mod p，其中r=2^p_bitcount
 * @param[in][out] a: 乘数，同时也是存放的输出结果
 * @param[in] b: 乘数
 * @param[in] p: 模数
 * @param[in] p_bitcount: 模数的比特数，确保是8的倍数，一般而言是2048或4096
 * @param[in] n_prime: 蒙哥马利需要的n'的值，和p长度相同
 */
void _mon_mod_mul_unit_128(char *a, char *b, char *p, int p_bitcount, char *n_prime) {
    // please first change '#define MON_MUL_BLOCK_LENGTH' as 128
    // MON_MUL_BLOCK_LENGTH = 128, p_bitcount = 2048, block_num = 16
    int block_num = (p_bitcount + MON_MUL_BLOCK_LENGTH - 1) / MON_MUL_BLOCK_LENGTH;
    int block_size = MON_MUL_BLOCK_LENGTH >> 3;
    int offset = block_size * block_num - BYTECOUNT(p_bitcount);
    int multi_block_size = block_size << 1;
    char *Q = new char[multi_block_size];
    char *temp = new char[multi_block_size];
    char *S = new char[block_size * (block_num + 2)];
    memset(S, 0x0, block_size * (block_num + 2));
    int carry_in_pre = 0;
    int carry_in_next = 0;
    bool carry_out = false;
    for (int i = 0; i <= block_num; i++) {
        // Q = (S0 + a0 * bi) * n_prime (mod 2^128)
        // Q = a0 * bi, a0 is the last block of a, bi is the (block_num - i)th block of b
        if (i < block_num) {
            _multi_128(a + (block_num - 1) * block_size - offset, b + (block_num - 1 - i) * block_size - offset, Q);
        }
        else {
            memset(Q, 0x0, multi_block_size);
        }
        // Q = Q mod 2^128
        memset(Q, 0x0, block_size);
        // temp = S0, S0 is the last block of S
        memset(temp, 0x0, block_size);
        memcpy(temp + block_size, S + (block_num + 1) * block_size, block_size);
        // Q = Q + S0
        _add_256(Q, temp, 0, Q, &carry_out);
        // Q = (Q mod 2^128) * n_prime
        _multi_128(Q + block_size, n_prime, Q);
        // Q = Q mod 2^128
        memset(Q, 0x0, block_size);
        // for j = 0, 1, ..., do: Sj = Sj + aj * bi + Q * pj
        carry_in_pre = 0;
        carry_in_next = 0;
        for (int j = 0; j < block_num; j++) {
            // temp = Q * pj, pj is the (block_num - j)th block of p
            _multi_128(Q + block_size, p + (block_num - 1 - j) * block_size, temp);
            // S = S + (temp << 128 * j)
            _add_256(temp, S + (block_num - j) * block_size, carry_in_pre, S + (block_num - j) * block_size, &carry_out);

            // process the carry
            carry_in_pre = carry_in_next;
            carry_in_next = (int)carry_out;
            // temp = aj * bi, aj is the (block_num - j)th block of a, bi is the (block_num - i)th block of b
            if (i < block_num) {
                _multi_128(a + (block_num - 1 - j) * block_size, b + (block_num - 1 - i) * block_size, temp);
            }
            else {
                memset(temp, 0x0, multi_block_size);
            }
            // S = S + (temp << 128 * j)
            _add_256(temp, S + (block_num - j) * block_size, 0, S + (block_num - j) * block_size, &carry_out);
            // process the carry
            carry_in_next = carry_in_next + (int)carry_out;
        }
        // process the last carry
        memset(temp, 0x0, multi_block_size);
        temp[block_size - 1] = (char)carry_in_next;
        temp[multi_block_size - 1] = (char)carry_in_pre;
        _add_256(temp, S, 0, S, &carry_out);

        // S = S / 2^128
        for (int j = block_size * (block_num + 2) - 1; j >= block_size; j--) {
            S[j] = S[j - block_size];
        }
        memset(S, 0x0, block_size);
    }

    // copy S to output
    for (int i = 0; i < block_size * block_num; i++) {
        a[i] = S[i + block_size * 2];
    }
    delete[] Q;
    delete[] temp;
    delete[] S;
}

/*
 * @brief: 蒙哥马利模乘单元，目标是输出 a*b*r^(-1) mod p，其中r=2^p_bitcount
 * @param[in][out] a: 乘数，同时也是存放的输出结果
 * @param[in] b: 乘数
 * @param[in] p: 模数
 * @param[in] p_bitcount: 模数的比特数，确保是8的倍数，一般而言是2048或4096
 * @param[in] n_prime: 蒙哥马利需要的n'的值，和p长度相同
 */
void _mon_mod_mul_unit_64(char *a, char *b, char *p, int p_bitcount, char *n_prime) {
    // MON_MUL_BLOCK_LENGTH = 64, p_bitcount = 2048, block_num = 32
    int block_num = (p_bitcount + MON_MUL_BLOCK_LENGTH - 1) / MON_MUL_BLOCK_LENGTH;
    int block_size = MON_MUL_BLOCK_LENGTH >> 3;
    int offset = block_size * block_num - BYTECOUNT(p_bitcount);
    int multi_block_size = block_size << 1;
    char *Q = new char[multi_block_size];
    char *temp = new char[multi_block_size];
    char *S = new char[block_size * (block_num + 2)];
    memset(S, 0x0, block_size * (block_num + 2));
    int carry_in_pre = 0;
    int carry_in_next = 0;
    bool carry_out = false;

    for (int i = 0; i <= block_num; i++) {
        // Q = (S0 + a0 * bi) * n_prime (mod 2^64)
        // Q = a0 * bi, a0 is the last block of a, bi is the (block_num - i)th block of b
        if (i < block_num) {
            _multi_64(a + (block_num - 1) * block_size - offset, b + (block_num - 1 - i) * block_size - offset, Q);
        }
        else {
            memset(Q, 0x0, multi_block_size);
        }

        // Q = Q mod 2^64
        memset(Q, 0x0, block_size);
        // temp = S0, S0 is the last block of S
        memset(temp, 0x0, block_size);
        memcpy(temp + block_size, S + (block_num + 1) * block_size, block_size);
        // Q = Q + S0
        _add_128(Q, temp, 0, Q, &carry_out);
        // Q = (Q mod 2^64) * n_prime
        _multi_64(Q + block_size, n_prime, Q);
        // Q = Q mod 2^64
        memset(Q, 0x0, block_size);

        // for j = 0, 1, ..., do: Sj = Sj + aj * bi + Q * pj
        carry_in_pre = 0;
        carry_in_next = 0;
        for (int j = 0; j < block_num; j++) {
            // temp = Q * pj, pj is the (block_num - j)th block of p
            _multi_64(Q + block_size, p + (block_num - 1 - j) * block_size, temp);
            // S = S + (temp << 64 * j)
            _add_128(temp, S + (block_num - j) * block_size, carry_in_pre, S + (block_num - j) * block_size, &carry_out);
            // process the carry
            carry_in_pre = carry_in_next;
            carry_in_next = (int)carry_out;

            // temp = aj * bi, aj is the (block_num - j)th block of a, bi is the (block_num - i)th block of b
            if (i < block_num) {
                _multi_64(a + (block_num - 1 - j) * block_size, b + (block_num - 1 - i) * block_size, temp);
            }
            else {
                memset(temp, 0x0, multi_block_size);
            }
            // S = S + (temp << 64 * j)
            _add_128(temp, S + (block_num - j) * block_size, 0, S + (block_num - j) * block_size, &carry_out);
            // process the carry
            carry_in_next = carry_in_next + (int)carry_out;
        }
        // process the last carry
        memset(temp, 0x0, multi_block_size);
        temp[block_size - 1] = (char)carry_in_next;
        temp[multi_block_size - 1] = (char)carry_in_pre;
        _add_128(temp, S, 0, S, &carry_out);

        // S = S / 2^64
        for (int j = block_size * (block_num + 2) - 1; j >= block_size; j--) {
            S[j] = S[j - block_size];
        }
        memset(S, 0x0, block_size);
    }

    // copy S to output
    for (int i = 0; i < block_size * block_num; i++) {
        a[i] = S[i + block_size * 2];
    }
    delete[] Q;
    delete[] temp;
    delete[] S;
}

/*
 * @brief: 蒙哥马利模幂单元，目标是输出 a^e*r mod p，其中r=2^p_bitcount
 * @param[in] a: 底数
 * @param[in] e: 指数
 * @param[in] e_bitcount: 指数的比特数，需要为8的倍数
 * @param[in] p: 模数
 * @param[in] p_bitcount: 模数的比特数，确保是8的倍数，一般而言是2048或4096
 * @param[in] n_prime: 蒙哥马利需要的n'的值，和p长度相同
 * @param[out] output: 输出的内容，硬件需要开辟一块单独的区域存放
 */
void _mon_mod_exp_unit(char *a, char *e, int e_bitcount, char *p, int p_bitcount, char *n_prime, char *r_square, char *output) {
    memset(output, 0, BYTECOUNT(p_bitcount));
    output[BYTECOUNT(p_bitcount) - 1] = 1;  // output存储的值是 2^p_bitcount
    _mon_mod_mul_unit_128(output, r_square, p, p_bitcount, n_prime);
    // 对e的每一位进行循环，由于C的操作单位是字节，所以拆成了两层循环
    for (int j = 0; j < BYTECOUNT(e_bitcount); ++j) {
        for (int k = 1; k <= BYTELEN; ++k) {
            _mon_mod_mul_unit_128(output, output, p, p_bitcount, n_prime);
            if (BIT(e[j], BYTELEN-k)) {
                _mon_mod_mul_unit_128(output, a, p, p_bitcount, n_prime);
            }
        }
    }
}

void mod_add(int count, char *a, char *b, char *p, int p_bitcount) {
    // 这count份计算是相互独立的，可并行进行
    for (int i = 0; i < count; ++i) {
        char *tmpa = OFFSET(a, p_bitcount, i);
        char *tmpb = OFFSET(b, p_bitcount, i);
        _mod_add_unit(tmpa, tmpb, p, p_bitcount);
    }
}

void mod_mul(int count, char *a, char *b, char *p, int p_bitcount, char *n_prime, char *r_square) {
    // one中存储的是1，由于每个操作都会用到，因此在这里生成一份给各操作共用
    char *one = new char[BYTECOUNT(p_bitcount)];
    memset(one, 0, BYTECOUNT(p_bitcount));
    one[BYTECOUNT(p_bitcount) - 1] = 1;
    // 这count份计算是相互独立的，可并行进行
    for (int i = 0; i < count; ++i) {
        char *tmpa = OFFSET(a, p_bitcount, i);
        char *tmpb = OFFSET(b, p_bitcount, i);
        _mon_mod_mul_unit_128(tmpa, r_square, p, p_bitcount, n_prime);
        _mon_mod_mul_unit_128(tmpb, r_square, p, p_bitcount, n_prime);
        _mon_mod_mul_unit_128(tmpa, tmpb, p, p_bitcount, n_prime);
        _mon_mod_mul_unit_128(tmpa, one, p, p_bitcount, n_prime);
    }
    delete [] one;
}

void mod_mul_ssl(int count, char *a, char *b, char *p, int p_bitcount, char *n_prime, char *r_square)
{
    BIGNUM *bn_a = BN_new();
    BIGNUM *bn_b = BN_new();
    BIGNUM *bn_p = BN_new();
    BIGNUM *bn_r = BN_new();
    BN_CTX * ctx = BN_CTX_new();
    BN_bin2bn((const unsigned char *)a, p_bitcount/8, bn_a);
    BN_bin2bn((const unsigned char *)b, p_bitcount/8, bn_b);
    BN_bin2bn((const unsigned char *)p, p_bitcount/8, bn_p);

    BN_mod_mul(bn_r, bn_a, bn_b, bn_p, ctx);
    BN_bn2binpad(bn_r, (unsigned char *)a, p_bitcount/8);

    BN_CTX_free(ctx);
    BN_free(bn_a);
    BN_free(bn_b);
    BN_free(bn_p);
    BN_free(bn_r);
}

void mod_mul_const(int count, char *a, char *b, char *p, int p_bitcount, char *n_prime, char *r_square) {
    // one中存储的是1，由于每个操作都会用到，因此在这里生成一份给各操作共用
    char *one = new char[BYTECOUNT(p_bitcount)];
    memset(one, 0, BYTECOUNT(p_bitcount));
    one[BYTECOUNT(p_bitcount) - 1] = 1;
    // 由于b是常数，因此可提前计算好b的蒙哥马利数并在后面复用
    _mon_mod_mul_unit_128(b, r_square, p, p_bitcount, n_prime);
    // 这count份计算是相互独立的，可并行进行
    for (int i = 0; i < count; ++i) {
        char *tmpa = OFFSET(a, p_bitcount, i);
        _mon_mod_mul_unit_128(tmpa, r_square, p, p_bitcount, n_prime);
        _mon_mod_mul_unit_128(tmpa, b, p, p_bitcount, n_prime);
        _mon_mod_mul_unit_128(tmpa, one, p, p_bitcount, n_prime);
    }
    delete [] one;
}

void mod_exp_ssl(int count, char *a, char *e, int e_bitcount, char *p, int p_bitcount, char *n_prime, char *r_square, char *output)
{
    BIGNUM *bn_a = BN_new();
    BIGNUM *bn_e = BN_new();
    BIGNUM *bn_p = BN_new();
    BIGNUM *bn_r = BN_new();
    BN_CTX * ctx = BN_CTX_new();
    BN_bin2bn((const unsigned char *)a, p_bitcount/8, bn_a);
    BN_bin2bn((const unsigned char *)e, e_bitcount/8, bn_e);
    BN_bin2bn((const unsigned char *)p, p_bitcount/8, bn_p);

    BN_mod_exp(bn_r, bn_a, bn_e, bn_p, ctx);
    BN_bn2binpad(bn_r, (unsigned char *)output, p_bitcount/8);

    BN_CTX_free(ctx);
    BN_free(bn_a);
    BN_free(bn_e);
    BN_free(bn_p);
    BN_free(bn_r);
}

void mod_exp(int count, char *a, char *e, int e_bitcount, char *p, int p_bitcount, char *n_prime, char *r_square, char *output) {
    // one中存储的是1，由于每个操作都会用到，因此在这里生成一份给各操作共用
    char *one = new char[BYTECOUNT(p_bitcount)];
    memset(one, 0, BYTECOUNT(p_bitcount));
    one[BYTECOUNT(p_bitcount) - 1] = 1;
    // 这count份计算是相互独立的，可并行进行
    for (int i = 0; i < count; ++i) {
        char *tmpa = OFFSET(a, p_bitcount, i);
        char *tmpe = OFFSET(e, e_bitcount, i);
        char *tmpo = OFFSET(output, p_bitcount, i);
        _mon_mod_mul_unit_128(tmpa, r_square, p, p_bitcount, n_prime);
        _mon_mod_exp_unit(tmpa, tmpe, e_bitcount, p, p_bitcount, n_prime, r_square, tmpo);
        _mon_mod_mul_unit_128(tmpo, one, p, p_bitcount, n_prime);
    }
    delete [] one;
}

void mod_const_exp(int count, char *a, char *e, int e_bitcount, char *p, int p_bitcount, char *n_prime, char *r_square, char *output) {
    // one中存储的是1，由于每个操作都会用到，因此在这里生成一份给各操作共用
    char *one = new char[BYTECOUNT(p_bitcount)];
    memset(one, 0, BYTECOUNT(p_bitcount));
    one[BYTECOUNT(p_bitcount) - 1] = 1;
    // 由于a是常数，因此可提前计算好a的蒙哥马利数并在后面复用
    _mon_mod_mul_unit_128(a, r_square, p, p_bitcount, n_prime);
    // 这count份计算是相互独立的，可并行进行
    for (int i = 0; i < count; ++i) {
        char *tmpe = OFFSET(e, e_bitcount, i);
        char *tmpo = OFFSET(output, p_bitcount, i);
        _mon_mod_exp_unit(a, tmpe, e_bitcount, p, p_bitcount, n_prime, r_square, tmpo);
        _mon_mod_mul_unit_128(tmpo, one, p, p_bitcount, n_prime);
    }
    delete [] one;
}

void mod_dot(int count, char *a, char *b, char *p, int p_bitcount, char *n_prime, char *r_square, char *output) {
    memset(output, 0, BYTECOUNT(p_bitcount));
    for (int i = 0; i < count; ++i) {
        char *tmpa = OFFSET(a, p_bitcount, i);
        char *tmpb = OFFSET(b, p_bitcount, i);
        _mon_mod_mul_unit_128(tmpb, r_square, p, p_bitcount, n_prime);
        _mon_mod_mul_unit_128(tmpa, tmpb, p, p_bitcount, n_prime);
        _mod_add_unit(output, tmpa, p, p_bitcount);
    }
}

void mod_exp_dot(int count, char *a, char *e, int e_bitcount, char *p, int p_bitcount, char *n_prime, char *r_square, char *output) {
    memset(output, 0, BYTECOUNT(p_bitcount));
    output[BYTECOUNT(p_bitcount) - 1] = 1;
    char tmpo[BYTECOUNT(p_bitcount)];
    for (int i = 0; i < count; ++i) {
        char *tmpa = OFFSET(a, p_bitcount, i);
        char *tmpe = OFFSET(e, e_bitcount, i);
        _mon_mod_mul_unit_128(tmpa, r_square, p, p_bitcount, n_prime);
        _mon_mod_exp_unit(tmpa, tmpe, e_bitcount, p, p_bitcount, n_prime, r_square, tmpo);
        _mon_mod_mul_unit_128(output, tmpo, p, p_bitcount, n_prime);
    }
}

/*
 * @brief: 密钥扩展
 * @param[in] key：初始16字节密钥
 * @param[out] exp_key：扩展176字节密钥，每32位已经倒置
 */
void aes_key_expand(const unsigned char *key, u32 *exp_key) {
    u32 temp;
    u32 *rk = exp_key;
    rk[0] = BUILD_UINT32(key[0], key[1], key[2], key[3]);
    rk[1] = BUILD_UINT32(key[4], key[5], key[6], key[7]);
    rk[2] = BUILD_UINT32(key[8], key[9], key[10], key[11]);
    rk[3] = BUILD_UINT32(key[12], key[13], key[14], key[15]);
    for (int i = 0; i < 10; i++) {
        temp  = rk[3];
        rk[4] = rk[0] ^
                (Te2[(temp >> 16) & 0xff] & 0xff000000) ^
                (Te3[(temp >>  8) & 0xff] & 0x00ff0000) ^
                (Te0[(temp) & 0xff] & 0x0000ff00) ^
                (Te1[(temp >> 24)] & 0x000000ff) ^
                rcon[i];
        rk[5] = rk[1] ^ rk[4];
        rk[6] = rk[2] ^ rk[5];
        rk[7] = rk[3] ^ rk[6];
        rk += 4;
    }
}

/*
 * @brief: 按128位为单元输出aes加密结果
 * @param[in]  in：输入128位明文数组
 * @param[in] exp_key：176字节轮密钥
 * @param[out] out：输出128位密文数组
 */
void aes_base_cipher(const unsigned char *in, unsigned char *out , const u32 *exp_key) {
    const u32 *rk = exp_key;
    unsigned char State[16];
    u32 s0, s1, s2, s3, t0, t1, t2, t3;
    s0 = BUILD_UINT32(in[0], in[1], in[2], in[3]) ^ rk[0];
    s1 = BUILD_UINT32(in[4], in[5], in[6], in[7]) ^ rk[1];
    s2 = BUILD_UINT32(in[8], in[9], in[10], in[11]) ^ rk[2];
    s3 = BUILD_UINT32(in[12], in[13], in[14], in[15]) ^ rk[3];
    for (int i= 1; i <= 5; i++ ) {
        t0 = Te0[s0 >> 24] ^ Te1[(s1 >> 16) & 0xff] ^ Te2[(s2 >>  8) & 0xff] ^ Te3[s3 & 0xff] ^ rk[8 * i - 4];
        t1 = Te0[s1 >> 24] ^ Te1[(s2 >> 16) & 0xff] ^ Te2[(s3 >>  8) & 0xff] ^ Te3[s0 & 0xff] ^ rk[8 * i - 3];
        t2 = Te0[s2 >> 24] ^ Te1[(s3 >> 16) & 0xff] ^ Te2[(s0 >>  8) & 0xff] ^ Te3[s1 & 0xff] ^ rk[8 * i - 2];
        t3 = Te0[s3 >> 24] ^ Te1[(s0 >> 16) & 0xff] ^ Te2[(s1 >>  8) & 0xff] ^ Te3[s2 & 0xff] ^ rk[8 * i - 1];
        if (i == 5) {
            break;
        }
        s0 = Te0[t0 >> 24] ^ Te1[(t1 >> 16) & 0xff] ^ Te2[(t2 >>  8) & 0xff] ^ Te3[t3 & 0xff] ^ rk[8 * i + 0];
        s1 = Te0[t1 >> 24] ^ Te1[(t2 >> 16) & 0xff] ^ Te2[(t3 >>  8) & 0xff] ^ Te3[t0 & 0xff] ^ rk[8 * i + 1];
        s2 = Te0[t2 >> 24] ^ Te1[(t3 >> 16) & 0xff] ^ Te2[(t0 >>  8) & 0xff] ^ Te3[t1 & 0xff] ^ rk[8 * i + 2];
        s3 = Te0[t3 >> 24] ^ Te1[(t0 >> 16) & 0xff] ^ Te2[(t1 >>  8) & 0xff] ^ Te3[t2 & 0xff] ^ rk[8 * i + 3];
    }
    s0 = (Te2[t0 >> 24] & 0xff000000) ^ (Te3[(t1 >> 16) & 0xff] & 0x00ff0000) ^
         (Te0[(t2 >> 8) & 0xff] & 0x0000ff00) ^  (Te1[t3 & 0xff] & 0x000000ff) ^ rk[40];
    s1 = (Te2[t1 >> 24] & 0xff000000) ^ (Te3[(t2 >> 16) & 0xff] & 0x00ff0000) ^
         (Te0[(t3 >> 8) & 0xff] & 0x0000ff00) ^ (Te1[t0 & 0xff] & 0x000000ff) ^ rk[41];
    s2 = (Te2[t2 >> 24] & 0xff000000) ^ (Te3[(t3 >> 16) & 0xff] & 0x00ff0000) ^
         (Te0[(t0 >> 8) & 0xff] & 0x0000ff00) ^ (Te1[t1 & 0xff] & 0x000000ff) ^ rk[42];
    s3 = (Te2[t3 >> 24] & 0xff000000) ^ (Te3[(t0 >> 16) & 0xff] & 0x00ff0000) ^
         (Te0[(t1 >> 8) & 0xff] & 0x0000ff00) ^ (Te1[t2 & 0xff] & 0x000000ff) ^ rk[43];
    PUTU32(out, s0);
    PUTU32(out + 4, s1);
    PUTU32(out + 8, s2);
    PUTU32(out + 12, s3);
}

/*
 * @brief: 按任意字节输出AES-CTR密文
 * @param[in]  iv：偏移量
 * @param[in]  exp_key：AES轮密钥
 * @param[in] length：所需密文字节长度
 * @param[out] out：输出密文数组
 */
void AES_CTR_encrypt_buffer(unsigned char *iv, u32  *exp_key, unsigned char* out, size_t length) {
    int runs = length / AES_BLOCK_SIZE;
    for (int i = 0; i < runs; ++i) {
        aes_base_cipher(iv, out + i * AES_BLOCK_SIZE, exp_key);
        for (int j = AES_BLOCK_SIZE - 1 ; j >= 0; j--) {
            if (++iv[j])
                break;
        }
    }
}

/*
 * @brief: 生成count个的随机数 r1, r2, ...
 * @param[in] count: 操作数的数目
 * @param[in] seed: 16字节
 * @param[in] bitcount: 每个随机数的比特位数
 * @param[in] output: 输出地址，有count个数，每个数是bitcount位
 */
void generate_pseudorandom(int count, unsigned char *seed, int bitcount, char *output) {
    int length = BYTECOUNT(count * bitcount);
    int temp_len = length;
    unsigned char *out = reinterpret_cast<unsigned char*>(output);
    unsigned char iv[16]={0};
    u32  exp_key[44];
    aes_key_expand(seed, exp_key);
    int padding_len = 0;
    if (length % AES_BLOCK_SIZE) {
        padding_len = AES_BLOCK_SIZE - (length % AES_BLOCK_SIZE);
    }
    bool need_padding = (padding_len != 0);
    int input_len1 = length;
    if (need_padding) {
        input_len1 = length - (AES_BLOCK_SIZE - padding_len);
    }
    if (length >= AES_BLOCK_SIZE) {
        AES_CTR_encrypt_buffer(iv, exp_key, out, input_len1);
    } else {
        input_len1 = 0;
    }
    if (need_padding) {
        unsigned char padding_buf[AES_BLOCK_SIZE] = {0};
        AES_CTR_encrypt_buffer(iv, exp_key, padding_buf, AES_BLOCK_SIZE);
        memcpy(out + input_len1, padding_buf, AES_BLOCK_SIZE - padding_len);
    }
}

/*
  * @brief: 多线程任务分配函数
  * @param[in] len：需要生成的随机数字节个数
  * @param[in]  working_thread_num：工作线程数
  * @param[in] [out] task：各线程分配字节数组
  * @param[in] [out] taskiv：各线程分配iv值
  */
void task_assign_thread(int len, int working_thread_num, int *task, unsigned char (*taskiv)[16]) {
    int temp = (len + AES_BLOCK_SIZE -1) / AES_BLOCK_SIZE;
    int temp1 = temp / working_thread_num;
    int  temp2 = temp % working_thread_num;
    int padding_len = 0;
    if (len % AES_BLOCK_SIZE) {
        padding_len = AES_BLOCK_SIZE - (len % AES_BLOCK_SIZE);
    }
    bool need_padding = (padding_len != 0);
    for (int i = 0; i < working_thread_num; i++){
        if (i < temp2){
            task[i] = (temp1 + 1) * AES_BLOCK_SIZE;
        }
        else {
            task[i] = temp1 * AES_BLOCK_SIZE;
        }
    }
    if (need_padding) {
        task[working_thread_num - 1] -= padding_len;
    }
    unsigned char ivtemp[16]={0};
    memcpy(taskiv[0], ivtemp, 16);
    int tempcount=0;
    for (int i = 1; i < working_thread_num; i++) {
        while (1) {
            for (int j = AES_BLOCK_SIZE - 1 ; j >= 0; j--) {
                if (++ivtemp[j])
                    break;
            }
            tempcount++;
            if (tempcount == task[i-1] / 16) {
                memcpy(taskiv[i], ivtemp, 16);
                tempcount = 0;
                break;
            }
        }
    }
}

/*
 * @brief: 多线程生成count个的随机数 r1, r2, ...
 * @param[in] count: 操作数的数目
 * @param[in] seed: 16字节种子
 * @param[in] bitcount: 每个随机数的比特位数
 * @param[in] thread_num: 可用线程数
 * @param[in] output: 输出地址，有count个数，每个数是bitcount位
 */
void generate_pseudorandom_multhread(int count, unsigned char *seed, int bitcount, char *output, int thread_num) {
    int length = BYTECOUNT(count * bitcount);
    unsigned char *out = reinterpret_cast<unsigned char*>(output);
    unsigned char iv[16]={0};
    u32  exp_key[44];
    aes_key_expand(seed, exp_key);
    /* 每个线程应至少分配到16字节的任务，否则不需要启动*/
    int working_thread_num = thread_num;
    if (AES_BLOCK_SIZE * thread_num > length) {
        working_thread_num = (length + AES_BLOCK_SIZE  - 1) / AES_BLOCK_SIZE;
    }
    /* 每个线程分配任务要保证生成顺序要与单线程一致，因此按16字节对齐*/
    int *task = new int[working_thread_num];
    unsigned char (*taskiv)[16] = new unsigned char[working_thread_num][16];
    task_assign_thread(length, working_thread_num, task, taskiv);
#pragma omp parallel for num_threads(working_thread_num)
    for (int i = 0; i < working_thread_num; i++) {
        size_t tasklen = task[i];
        int threadout_offset = 0;
        for (int j =0 ; j < i; j++){
            threadout_offset += task[j];
        }
        if (i < working_thread_num  - 1) {
            AES_CTR_encrypt_buffer(taskiv[i], exp_key, out + threadout_offset, tasklen);
        } else {
            int padding_len = 0;
            if (length % AES_BLOCK_SIZE) {
                padding_len = AES_BLOCK_SIZE - (length % AES_BLOCK_SIZE);
            }
            bool need_padding = (padding_len != 0);
            int input_len1 = tasklen;
            if (need_padding) {
                input_len1 = tasklen - (AES_BLOCK_SIZE - padding_len);
            }
            if (tasklen >= AES_BLOCK_SIZE) {
                AES_CTR_encrypt_buffer(taskiv[i], exp_key, out + threadout_offset, input_len1);
            } else {
                input_len1 = 0;
            }
            if (need_padding) {
                unsigned char padding_buf[AES_BLOCK_SIZE] = {0};
                AES_CTR_encrypt_buffer(iv, exp_key, padding_buf, AES_BLOCK_SIZE);
                memcpy(out + threadout_offset + input_len1, padding_buf, AES_BLOCK_SIZE - padding_len);
            }
        }
    }
    delete [] task;
    delete [] taskiv;
}

/*
 * @brief: 大整数右移
 * @param[in] a：输入大整数
 * @param[in] operating_bytes：处理字节大小
 */
void right_shift_1bit (char *a, int operating_bytes) {
    unsigned char* b = (unsigned char*) a;
    for (int i = operating_bytes -1; i > 0; i--) {
        b[i] = b[i] >> 1;
        if (b[i - 1] & LSB)
            b[i] = b[i] | MSB;
    }
    if (b[0] & MSB)
        b[0] = (b[0] >> 1) | MSB;
    else
        b[0] = b[0] >> 1;
}

/*
 * @brief: 判断大整数是否为偶数
 * @param[in] a：输入大整数
 * @param[in] operating_bytes：处理字节大小
 * @param[out]：若a为偶数返回true，否则false
 */
bool judge_bignum_even(char *a, int operating_bytes) {
    if (a[operating_bytes - 1] & LSB)
        return false;
    return true;
}

/*
 * @brief: 判断大整数是否为1
 * @param[in] a：输入大整数
 * @param[in] operating_bytes：处理字节大小
 * @param[out]：若a为1返回true，否则false
 */
bool judge_bignum_equal_one(char *a, int operating_bytes) {
    if (a[operating_bytes - 1] != 1) {
        return false;
    } else {
        for (int i = 0; i < operating_bytes -1; i++) {
            if (a[i] != 0)
                return false;
        }
        return true;
    }
}

/*
 * @brief: 判断大整数是否为0
 * @param[in] a：输入大整数
 * @param[in] operating_bytes：处理字节大小
 * @param[out]：若a为0返回true，否则false
 */
bool judge_bignum_equal0(char *a, int operating_bytes) {
    for (int i = 0; i < operating_bytes ; i++) {
        if (a[i] != 0)
            return false;
    }
    return true;
}

/*
 * @brief: 判断大整数a,b是否满足a > b
 * @param[in] a：输入大整数
 * @param[in] b：输入大整数
 * @param[in] operating_bytes：处理字节大小
 * @param[out]：若a>b返回true，否则false
 */
bool judge_bignum_larger(char *a, char *b, int operating_bytes) {
    bool a_lager_0 = (a[0] >= 0);
    bool b_lager_0 = (b[0] >= 0);
    if (a_lager_0  && !b_lager_0) // if a >= 0, b < 0
        return true;
    else if (!a_lager_0 && b_lager_0) // if a < 0, b >= 0
        return false;
    else if (a_lager_0 && b_lager_0) { // if a >= 0, b >= 0
        unsigned char temp_achar;
        unsigned char temp_bchar;
        for (int i = 0; i <operating_bytes ; i++) {
            temp_achar = a[i];
            temp_bchar = b[i];
            if (temp_achar > temp_bchar) {
                return true;
            }
            if (temp_achar < temp_bchar) {
                return false;
            }
        }
        return false;
    } else { // if a < 0 b < 0
        if (a[0] > b[0])
            return true;
        unsigned char temp_achar;
        unsigned char temp_bchar;
        for (int i = 1; i <operating_bytes; i++) {
            temp_achar = a[i];
            temp_bchar = b[i];
            if (temp_achar < temp_bchar) {
                return true;
            }
            if (temp_achar > temp_bchar) {
                return false;
            }
        }
        return false;
    }
}

/*
 * @brief: 大整数相加，非模加
 * @param[in] [out] a：输入大整数，同时存放相加结果
 * @param[in] b：输入大整数
 * @param[in] operating_bytes：处理字节大小
 */
void bignum_add(char *a, char *b, int operating_bytes) {
    if (operating_bytes <= 16) {
        __uint128_t a_temp = 0;
        __uint128_t b_temp = 0;
        _decode_char_to_int(a, operating_bytes, &a_temp);
        _decode_char_to_int(b, operating_bytes, &b_temp);
        a_temp += b_temp;
        _encode_int_to_char(&a_temp, operating_bytes, a);
        return;
    }
    int block_size = 16;
    int block_count = (operating_bytes + block_size - 1) / block_size;
    // 在加减法实现中，carry_in和carry_out不能指向同一块空间
    int carry_in = 0;
    bool carry_out = false;
    for (int i = block_count - 1; i >= 0; i--) {
        _add_128(a + i * block_size, b + i * block_size, carry_in,
                 a + i * block_size, &carry_out);
        // process the carry
        carry_in = (int)carry_out;
    }
}

/*
 * @brief: 大整数取相反数
 * @param[in] a：输入大整数
 * @param[in] operating_bytes：处理字节大小
 * @param[out] a_neg：a的相反数结果
 */
void bignum_negtive(char *a, char *a_neg, int operating_bytes) {
    int block_size = 16;
    int block_count = (operating_bytes + block_size - 1) / block_size;
    __uint128_t a_temp = 0;
    for (int i = 0; i < block_count; i++) {
        _decode_char_to_int(a + i * block_size , sizeof(a_temp), &a_temp);
        a_temp = ~a_temp;
        _encode_int_to_char(&a_temp, sizeof(a_temp), a_neg + i * block_size);
    }
    char *temp1 = new char[operating_bytes];
    memset(temp1, 0, operating_bytes);
    temp1[operating_bytes - 1] = 1;
    bignum_add(a_neg, temp1, operating_bytes);
    delete [] temp1;
}

/*
 * @brief: 求解模逆单元，目标是输出 a^{-1} mod p
 * @param[in] [out] a:  求逆，也是输出结果，若不存在模逆输出结果0
 * @param[in] p: 模数
 * @param[in] p_bitcount: 模数的比特数，确保是64的倍数，一般而言是2048或4096
 */
void _mod_inverse_unit(char *a, char *p, int p_bitcount) {
    int operating_bytes = BYTECOUNT(p_bitcount);
    if (judge_bignum_even(a, operating_bytes) && judge_bignum_even(p, operating_bytes)) {
        memset(a, 0, operating_bytes);
        std::cout<<"a and p are both even, no inverse"<<std::endl;
        return;
    }
    char *x1 = new char[operating_bytes];
    memset(x1, 0, operating_bytes);
    x1[operating_bytes - 1] = 1;

    char *x2 = new char[operating_bytes];
    memset(x2, 0, operating_bytes);

    char *x3 = new char[operating_bytes];
    memcpy(x3, a, operating_bytes);

    char *y1 = new char[operating_bytes];
    memset(y1, 0, operating_bytes);

    char *y2 = new char[operating_bytes];
    memset(y2, 0, operating_bytes);
    y2[operating_bytes -1] = 1;

    char *y3 = new char[operating_bytes];
    memcpy(y3, p, operating_bytes);

    char *a_neg =new char[operating_bytes];
    bignum_negtive(a, a_neg, operating_bytes);

    char *p_neg =new char[operating_bytes];
    bignum_negtive(p, p_neg, operating_bytes);

    while ( !judge_bignum_equal_one(x3, operating_bytes) && !judge_bignum_equal_one(y3, operating_bytes)){
        if (judge_bignum_equal0(x3, operating_bytes) || judge_bignum_equal0(y3, operating_bytes)){
            memset(a, 0, operating_bytes);
            std::cout<<"gcd != 1, no inverse"<<std::endl;
            return;
        }

        if (judge_bignum_even(x3, operating_bytes)) { // if x3 is even
            if (judge_bignum_even(x1, operating_bytes) && judge_bignum_even(x2, operating_bytes)) {
                right_shift_1bit(x1, operating_bytes);
                right_shift_1bit(x2, operating_bytes);
                right_shift_1bit(x3, operating_bytes);
            } else {
                bignum_add(x1, p, operating_bytes);
                bignum_add(x2, a_neg, operating_bytes);
            }
        } else if (judge_bignum_even(y3, operating_bytes)) { // if y3 is even
            if (judge_bignum_even(y1, operating_bytes) && judge_bignum_even(y2, operating_bytes)) {
                right_shift_1bit(y1, operating_bytes);
                right_shift_1bit(y2, operating_bytes);
                right_shift_1bit(y3, operating_bytes);
            } else {
                bignum_add(y1, p, operating_bytes);
                bignum_add(y2, a_neg, operating_bytes);
            }
        } else if (judge_bignum_larger(x3, y3, operating_bytes)) { // if both x3, y3 is not even, and x3 > y3
            char *neg_temp1 = new char[operating_bytes];
            bignum_negtive(y1, neg_temp1, operating_bytes);
            char *neg_temp2 = new char[operating_bytes];
            bignum_negtive(y2, neg_temp2, operating_bytes);
            char *neg_temp3 = new char[operating_bytes];
            bignum_negtive(y3, neg_temp3, operating_bytes);

            bignum_add(x1, neg_temp1, operating_bytes);
            bignum_add(x2, neg_temp2, operating_bytes);
            bignum_add(x3, neg_temp3, operating_bytes);
            delete [] neg_temp1;
            delete [] neg_temp2;
            delete [] neg_temp3;
        } else {
            char *neg_temp1 = new char[operating_bytes];
            bignum_negtive(x1, neg_temp1, operating_bytes);
            char *neg_temp2 = new char[operating_bytes];
            bignum_negtive(x2, neg_temp2, operating_bytes);
            char *neg_temp3 = new char[operating_bytes];
            bignum_negtive(x3, neg_temp3, operating_bytes);

            bignum_add(y1, neg_temp1, operating_bytes);
            bignum_add(y2, neg_temp2, operating_bytes);
            bignum_add(y3, neg_temp3, operating_bytes);
            delete [] neg_temp1;
            delete [] neg_temp2;
            delete [] neg_temp3;
        }
    } // end while

    if (judge_bignum_equal_one(x3, operating_bytes)) {
        memcpy(a, x1, operating_bytes);
    } else {
        memcpy(a, y1, operating_bytes);
    }
    while (a[0] < 0) {  //  if a < 0, output += p
        bignum_add(a, p, operating_bytes);
    }
    while (judge_bignum_larger(a,p,operating_bytes)) {  //  if a > p, output -= p
        bignum_add(a, p_neg, operating_bytes);
    }
    delete [] x1;
    delete [] x2;
    delete [] x3;
    delete [] y1;
    delete [] y2;
    delete [] y3;
    delete [] a_neg;
    delete [] p_neg;
    return;
}

/*
 * @brief: 实现模数上求逆，目标是输出 a1^{-1} mod p, a2^{-1} mod p, ...
 * @param[in] count: 操作数的数目
 * @param[in][out] a: 包含count个值，每个有p_bitcount个比特，也是输出存放的位置
 * @param[in] p: 模数
 * @param[in] p_bitcount: 模数的比特数，确保是64的倍数，一般而言是2048或4096
 */
void mod_inverse(int count, char *a, char *p, int p_bitcount) {
    // 这count份计算是相互独立的，可并行进行
    for (int i = 0; i < count; i++) {
        char *tmpa = OFFSET(a, p_bitcount, i);
        _mod_inverse_unit(tmpa, p, p_bitcount);
    }
}

void operation(char opcode, int count, int p_bitcount, char *input, char **output) {
    int offset = 0;
    char *p, *n_prime, *r_square;

    if (FOURBITS(opcode) <= 8) {  // 所有已支持的接口都需要获取模数p的值
        // 读出模数p的值
        p = input + offset; offset += BYTECOUNT(p_bitcount);
    }
    if (1 <= FOURBITS(opcode) && FOURBITS(opcode) <= 7) {  // 中间某些操作需要这些值
        // 读出蒙哥马利计算时需要的n'的值
        n_prime = input + offset; offset += BYTECOUNT(MON_MUL_BLOCK_LENGTH);
        // 读出蒙哥马利计算时需要的r^2的值
        r_square = input + offset; offset += BYTECOUNT(p_bitcount);
    }

    if (FOURBITS(opcode) == 0) {  // 操作码为0，代表向量模加
        char *a = input + offset; offset += BYTECOUNT(p_bitcount) * count;
        char *b = input + offset; offset += BYTECOUNT(p_bitcount) * count;
        mod_add(count, a, b, p, p_bitcount);
        // 模加的输出复用了a的空间
        *output = a;
    } else if (FOURBITS(opcode) == 1) {  // 操作码为1，代表向量模乘
        char *a = input + offset; offset += BYTECOUNT(p_bitcount) * count;
        char *b = input + offset; offset += BYTECOUNT(p_bitcount) * count;
        mod_mul(count, a, b, p, p_bitcount, n_prime, r_square);
        // 模乘的输出复用了a的空间
        *output = a;
    } else if (FOURBITS(opcode) == 2) {  // 操作码为2，代表向量模乘，不过第二个乘数是个常量而非向量
        char *a = input + offset; offset += BYTECOUNT(p_bitcount) * count;
        char *b = input + offset; offset += BYTECOUNT(p_bitcount);
        mod_mul_const(count, a, b, p, p_bitcount, n_prime, r_square);
        // 模乘的输出复用了a的空间
        *output = a;
    } else if (FOURBITS(opcode) ==3) {  // 操作码为3，代表向量模幂
        char *a = input + offset; offset += BYTECOUNT(p_bitcount) * count;
        // 读出指数的位长度，一般而言不会超过4096，所以两个字节（即16个比特）足够
        short int e_bitcount = *(reinterpret_cast<short int*>(input + offset)); offset += sizeof(short int);
        char *e = input + offset; offset += BYTECOUNT(e_bitcount) * count;
        // 模幂的输出需要在硬件中额外开辟空间
        *output = new char[BYTECOUNT(p_bitcount) * count];
        mod_exp(count, a, e, e_bitcount, p, p_bitcount, n_prime, r_square, *output);
    } else if (FOURBITS(opcode) == 4) {  // 操作码为4，代表向量模幂，不过底数是个常量而非向量
        char *a = input + offset; offset += BYTECOUNT(p_bitcount);
        // 读出指数的位长度，一般而言不会超过4096，所以两个字节（即16个比特）足够
        short int e_bitcount = *(reinterpret_cast<short int*>(input + offset)); offset += sizeof(short int);
        char *e = input + offset; offset += BYTECOUNT(e_bitcount) * count;
        // 模幂的输出需要在硬件中额外开辟空间
        *output = new char[BYTECOUNT(p_bitcount) * count];
        mod_const_exp(count, a, e, e_bitcount, p, p_bitcount, n_prime, r_square, *output);
    } else if (FOURBITS(opcode) == 5) {  // 操作码为5，代表同态乘法点积
        char *a = input + offset; offset += BYTECOUNT(p_bitcount) * count;
        char *b = input + offset; offset += BYTECOUNT(p_bitcount) * count;
        // 输出需要在硬件中额外开辟空间
        *output = new char[BYTECOUNT(p_bitcount)];
        mod_dot(count, a, b, p, p_bitcount, n_prime, r_square, *output);
    } else if (FOURBITS(opcode) == 6) {  // 操作码为6，代表同态指数点积
        char *a = input + offset; offset += BYTECOUNT(p_bitcount) * count;
        // 读出指数的位长度，一般而言不会超过4096，所以两个字节（即16个比特）足够
        short int e_bitcount = *(reinterpret_cast<short int*>(input + offset)); offset += sizeof(short int);
        char *e = input + offset; offset += BYTECOUNT(e_bitcount) * count;
        // 输出需要在硬件中额外开辟空间
        *output = new char[BYTECOUNT(p_bitcount)];
        mod_exp_dot(count, a, e, e_bitcount, p, p_bitcount, n_prime, r_square, *output);
    } else if (FOURBITS(opcode) == 7) {  // 操作码为7，代表批量求逆
        // 保留操作码，待后续实现
    } else if (FOURBITS(opcode) == 8) {  // 操作码为8，代表批量生成随机数
        // 保留操作码，待后续实现
    } else {
        std::cout << "wrong operation code " << int(FOURBITS(opcode)) << std::endl;
    }
}
